#include "TFile.h"
#include "TH2F.h"
#include "TGraph2D.h"
#include "TGraphDelaunay.h"
#include "TCanvas.h"
#include <iostream>
void createPDF()
{
    const unsigned int numFiles = 9;
    const unsigned int MAXPOINTS = 10000;
    const unsigned int newBins = 200;
    TString str_atoms[numFiles] = {"Ar", "C", "F", "Ge", "He", "Kr", "Ne", "Si", "Xe"};
    TString str_files[numFiles] = {"Ar+.root", "C+.root", "F+.root", "Ge+.root", "He+.root", "Kr+.root", "Ne+.root", "Si+.root", "Xe+.root"};

    TFile *newFile = new TFile("MigdalEnergyPDF.root", "RECREATE");

    for (unsigned int ia = 0; ia < numFiles; ++ia)
    {

        TFile *file = new TFile(str_files[ia]);
        TH2F *hist = (TH2F *)file->Get("PDF");
        TH1 *hist_x = hist->ProjectionX("PDF_x");
        TH1 *hist_y = hist->ProjectionY("PDF_y");
        unsigned int numPoints = 0;
        unsigned int nBinX = hist->GetNbinsX();
        unsigned int nBinY = hist->GetNbinsY();

        Double_t x_array[MAXPOINTS] = {0.};
        Double_t y_array[MAXPOINTS] = {0.};
        Double_t z_array[MAXPOINTS] = {0.};
        for (unsigned int ix = 1; ix < nBinX + 1; ++ix)
        {
            for (unsigned int iy = 1; iy < nBinY + 1; ++iy)
            {
                Double_t val = hist->GetBinContent(ix, iy);
                if (val > 0)
                {
                    x_array[numPoints] = hist_x->GetBinCenter(ix);
                    y_array[numPoints] = hist_y->GetBinCenter(iy);
                    z_array[numPoints] = val;
                    ++numPoints;
                }
            }
        }
        TGraph2D *graph = new TGraph2D(numPoints, x_array, y_array, z_array);
        TGraphDelaunay graphDe(graph);
        Double_t xmin = graph->GetXmin();
        Double_t xmax = graph->GetXmax();
        Double_t ymin = graph->GetYmin();
        Double_t ymax = graph->GetYmax();
        TH2F *newHist = new TH2F(str_atoms[ia], "", newBins, xmin, xmax, newBins, ymin, ymax);
        TH1 *newHist_x = newHist->ProjectionX("new_x");
        TH1 *newHist_y = newHist->ProjectionY("new_y");
        for (unsigned int ix = 1; ix < newBins + 1; ++ix)
        {
            for (unsigned int iy = 1; iy < newBins + 1; ++iy)
            {
                Double_t x = newHist_x->GetBinCenter(ix);
                Double_t y = newHist_y->GetBinCenter(iy);
                Double_t z = graphDe.ComputeZ(x, y);
                newHist->SetBinContent(ix, iy, z);
            }
        }
        newFile->cd();
        newHist->Write();
        numPoints = 0;
        file->Close();
        std::cout << "PDF of atom " << str_atoms[ia] << " saved" << std::endl;
    }
    newFile->Close();
}